# Load Necessary Packages
library(glmnet)
library(foreign)
library(estimatr)
library(dplyr)
library(tidyr)
library(bcf)
library(dbarts)
library(texreg)
library(Hmisc)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(car)
library(plotrix)
library(ggpubr)
library(effectsize)
library(MASS)
library(rms)
library(pwr)
library(DeclareDesign)
library(estimatr)
library(haven)
library(TOSTER)

# Load Data
setwd("/users/josephphillips/Dropbox/COVID-19/Data/")
data <- read.dta("US_pooled_complete_v12.dta")

# Load glmnet function
run_LASSO <- function(data, yvar){
  
  set.seed(2334)
  
  var_quo <- enquo(yvar)
  
  ## select x variables 
  dat <- data %>%
    dplyr::select(college, agegroup, male, married, frequent_church, region, pid3, ideo7, high_incidence,crt, knowledge,nonwhite, polinterest,health_trust,media_trust)
  
  ## select y-variable
  out1 <- data %>% 
    dplyr::select(!!var_quo) %>% pull() # save as a vector instead of a data-frame to use model.matrix()
  
  merged1 <- cbind(out1, dat) %>% na.omit() ## merge and delete NA
  x1 <- model.matrix(out1 ~ ., data = merged1) ## make model matrix
  y1 <- merged1$out1 ## select y-variable from NA-deleted dataset
  mod1 <- cv.glmnet(x1, y1, alpha = 1, family = "gaussian") ## run lasso model
  coef <- coef(mod1) ## get coef  
  
  list(merged1, x1, y1, mod1, coef)
  
}


# Clean Controls
data$male <- as.numeric(data$male) -1
data$married <- as.numeric(data$married) - 1
data$frequent_church <- as.numeric(data$frequent_church) - 1
data$pid3 <- ifelse(data$democrat==1,0,ifelse(data$republican==1,2,ifelse(data$pid7=="Independent",1,NA)))
data$ideo7 <- as.numeric(data$ideo7)
data$high_incidence <- as.numeric(data$high_incidence) - 1
data$lives_highincidence <- as.numeric(data$lives_highincidence) - 1
data$nonwhite <- as.numeric(data$nonwhite) - 1
data$polinterest <- as.numeric(data$polinterest)
data$factcheck_treatment <- as.numeric(data$factcheck_treatment) - 1
data$W3factcheck_treatment <- as.numeric(data$W3factcheck_treatment) -1
data$approve_trmp <- as.numeric(data$approve_trmp) - 1
data$conspiracy <- data$conspiracy-1
data$health_trust <- data$health_trust-1
data$media_trust <- data$media_trust-1
data$therm_trump <- data$therm_trump/100

# Clean DVs
data$wear_mask <- 6-(data$wear_mask) ## It was reverse-coded
data$mask_effective <- as.numeric(data$mask_effective)
## Note: with affective polarization, the party affect variables are actually mislabeled
data$affective_polarization <- ifelse(data$pid3==0,data$feeling_Rep-data$feeling_Dem,ifelse(data$pid3==2,data$feeling_Dem-data$feeling_Rep,NA))

# Clean Treatment
data$american_treatment <- ifelse(data$exposure_treat=="Americans",1,0)
data$copartisan_treatment <- ifelse((data$pid3==0 & data$exposure_treat=="Democrats") | (data$pid3==2 & data$exposure_treat=="Republicans"),1,0)
data$contrapartisan_treatment <- ifelse((data$pid3==0 & data$exposure_treat=="Republicans") | (data$pid3==2 & data$exposure_treat=="Democrats"),1,0)

# Create underestimation/overestimation in mask wearing
data$mask_allthetime <- ifelse(data$W2_covid19_actions_2=="Most of the time" | data$W2_covid19_actions_2=="All of the time",1,0)
data$underestimates_masks_american <- ifelse(data$wear_mask_pub<70,1,0)
data$overestimates_masks_american <- ifelse(data$wear_mask_pub>=90,1,0)
data$underestimates_masks_copartisans <- ifelse(data$pid3==2 & data$wear_mask_Rep<47,1,ifelse(data$pid3==0 & data$wear_mask_Dem<80,1,ifelse(data$pid3==1,NA,0)))
data$overestimates_masks_copartisans <- ifelse(data$pid3==2 & data$wear_mask_Rep>56,1,ifelse(data$pid3==0 & data$wear_mask_Dem>98,1,ifelse(data$pid3==1,NA,0)))
data$underestimates_masks_contrapartisans <- ifelse(data$pid3==0 & data$wear_mask_Rep<47,1,ifelse(data$pid3==2 & data$wear_mask_Dem<80,1,ifelse(data$pid3==1,NA,0)))
data$overestimates_masks_contrapartisans <- ifelse(data$pid3==0 & data$wear_mask_Rep>56,1,ifelse(data$pid3==2 & data$wear_mask_Dem>98,1,ifelse(data$pid3==1,NA,0))) ## Note: It's impossible to overestimate Democrat mask wearing by 10%
data$underestimates_masks_republicans <- ifelse(data$wear_mask_Rep<47,1,0)
data$overestimates_masks_republicans <- ifelse(data$wear_mask_Rep>56,1,0)
data$underestimates_masks_democrats <- ifelse(data$wear_mask_Dem<80,1,0)
data$overestimates_masks_democrats <- ifelse(data$wear_mask_Dem>98,1,0)

# Lasso Regression
df <- data
wear_lasso <- run_LASSO(df,wear_mask) ## male, pid3, ideo7, health_trust
effective_lasso <- run_LASSO(df,mask_effective) ## pid3, ideo7, health_trust, media_trust
polariz_lasso <- run_LASSO(df,affective_polarization) ## ideo7, polinterest

# Table 1
## Report Regularly Wearing Masks
weighted.mean(data$mask_allthetime,weight=weight_genpop_pulse_high_inciden,na.rm=T) ## 79.54% of Americans
weighted.mean(subset(data,pid3==0)$mask_allthetime,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T) ## 93.14% of Democrats
weighted.mean(subset(data,pid3==2)$mask_allthetime,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T) ## 60.11% of Republicans
## Estimates of American Mask-Wearing
weighted.mean(data$wear_mask_pub,data$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==0)$wear_mask_pub,subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==2)$wear_mask_pub,subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)
## Estimates of Democrat Mask-Wearing
weighted.mean(data$wear_mask_Dem,data$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==0)$wear_mask_Dem,subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==2)$wear_mask_Dem,subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)
## Estimates of Republican Mask-Wearing
weighted.mean(data$wear_mask_Rep,data$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==0)$wear_mask_Rep,subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)
weighted.mean(subset(data,pid3==2)$wear_mask_Rep,subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)
## Underestimates American Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$underestimates_masks_american,weight=data$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$underestimates_masks_american,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$underestimates_masks_american,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Overestimates American Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$overestimates_masks_american,weight=data$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$overestimates_masks_american,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$overestimates_masks_american,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Underestimates Democratic Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$underestimates_masks_democrats,weight=weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$underestimates_masks_democrats,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$underestimates_masks_democrats,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Overestimates Democratic Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$overestimates_masks_democrats,weight=weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$overestimates_masks_democrats,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$overestimates_masks_democrats,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Underestimates Republican Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$underestimates_masks_republicans,weight=weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$underestimates_masks_republicans,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$underestimates_masks_republicans,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Overestimates Republican Mask-Wearing by 10% or More (CHECK)
weighted.mean(data$overestimates_masks_republicans,weight=weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==0)$overestimates_masks_republicans,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)*100
weighted.mean(subset(data,pid3==2)$overestimates_masks_republicans,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)*100
## Mask-Wearing Intentons
weighted.mean(data$wear_mask,weight=data$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(data$wear_mask,weight=data$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==0)$wear_mask,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(subset(data,pid3==0)$wear_mask,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==2)$wear_mask,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(subset(data,pid3==2)$wear_mask,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden))
## Perceived Mask Effectiveness
weighted.mean(data$mask_effective,weight=data$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(data$mask_effective,weight=data$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==0)$mask_effective,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(subset(data,pid3==0)$mask_effective,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==2)$mask_effective,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden)
sqrt(wtd.var(subset(data,pid3==2)$mask_effective,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden))
## Affective Polarization
weighted.mean(data$affective_polarization,weight=data$weight_genpop_pulse_high_inciden,na.rm=T)
sqrt(wtd.var(data$affective_polarization,weight=data$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==0)$affective_polarization,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden,na.rm=T)
sqrt(wtd.var(subset(data,pid3==0)$affective_polarization,weight=subset(data,pid3==0)$weight_genpop_pulse_high_inciden))
weighted.mean(subset(data,pid3==2)$affective_polarization,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden,na.rm=T)
sqrt(wtd.var(subset(data,pid3==2)$affective_polarization,weight=subset(data,pid3==2)$weight_genpop_pulse_high_inciden))

# Figure 1
## Regression
Reg.h1a.partisans <- lm_robust(wear_mask~american_treatment+copartisan_treatment+contrapartisan_treatment+male+pid3+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2)
## Standardized Effect Size for American Treatment
Reg.h1a.partisans$coefficients[2]/sd(data$wear_mask,na.rm=T)
## Equivalence Tests
data$treatment <- as.factor(ifelse(data$american_treatment==1,"American",ifelse(data$copartisan_treatment==1,"Co-Partisan",ifelse(data$contrapartisan_treatment==1,"Out-Partisan","Control"))))
data.partisans <- subset(data,pid3==0 | pid3==2)
diffmeans.h1a.partisans <- data.partisans %>% group_by(treatment) %>% summarise(mean=mean(wear_mask,na.rm=T), sd=sd(wear_mask,na.rm=T), n=n())
### Control vs. Co-Partisan
TOSTtwo.raw(m1=diffmeans.h1a.partisans$mean[2],m2=diffmeans.h1a.partisans$mean[3],sd1=diffmeans.h1a.partisans$sd[2],sd2=diffmeans.h1a.partisans$sd[3],n1=diffmeans.h1a.partisans$n[2],n2=diffmeans.h1a.partisans$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F)
### Control vs. Out-Partisan
TOSTtwo.raw(m1=diffmeans.h1a.partisans$mean[2],m2=diffmeans.h1a.partisans$mean[4],sd1=diffmeans.h1a.partisans$sd[2],sd2=diffmeans.h1a.partisans$sd[4],n1=diffmeans.h1a.partisans$n[2],n2=diffmeans.h1a.partisans$n[4],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Out-Partisan vs. American
TOSTtwo.raw(m1=diffmeans.h1a.partisans$mean[4],m2=diffmeans.h1a.partisans$mean[1],sd1=diffmeans.h1a.partisans$sd[4],sd2=diffmeans.h1a.partisans$sd[1],n1=diffmeans.h1a.partisans$n[4],n2=diffmeans.h1a.partisans$n[1],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Out-Patisan vs. Co-Partisan
TOSTtwo.raw(m1=diffmeans.h1a.partisans$mean[4],m2=diffmeans.h1a.partisans$mean[3],sd1=diffmeans.h1a.partisans$sd[4],sd2=diffmeans.h1a.partisans$sd[3],n1=diffmeans.h1a.partisans$n[4],n2=diffmeans.h1a.partisans$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
## Comparing Effect Sizes
linearHypothesis(Reg.h1a.partisans,"american_treatment=copartisan_treatment")
linearHypothesis(Reg.h1a.partisans,"american_treatment=contrapartisan_treatment")
linearHypothesis(Reg.h1a.partisans,"copartisan_treatment=contrapartisan_treatment")
## Create Figure
Reg.h1a.partisans.strict <- confint(Reg.h1a.partisans,level=.995)
newdata1 <- data.frame(treatment=factor(c("American","Copartisan","Contrapartisan"),levels=c("Contrapartisan","Copartisan","American")),
                       effect=c(Reg.h1a.partisans$coefficients[2],Reg.h1a.partisans$coefficients[3],Reg.h1a.partisans$coefficients[4]),
                       lci=c(Reg.h1a.partisans$conf.low[2],Reg.h1a.partisans$conf.low[3],Reg.h1a.partisans$conf.low[4]),
                       uci=c(Reg.h1a.partisans$conf.high[2],Reg.h1a.partisans$conf.high[3],Reg.h1a.partisans$conf.high[4]),
                       lci_strict=c(Reg.h1a.partisans.strict[2,1],Reg.h1a.partisans.strict[3,1],Reg.h1a.partisans.strict[4,1]),
                       uci_strict=c(Reg.h1a.partisans.strict[2,2],Reg.h1a.partisans.strict[3,2],Reg.h1a.partisans.strict[4,2]))

p.h1a <- ggplot(newdata1,aes(x=treatment,y=effect)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=.8,width=0) + geom_errorbar(aes(ymin=lci_strict,ymax=uci_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53))

# Figure 2
## Regression
Reg.h1b.partisans <- lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans)+male+pid3+ideo7+health_trust,data=data)
summary(margins(Reg.h1b.partisans,at=list(underestimates_masks_american=1))) ## Marginal effect for underestimators
## Figure
h1b.underestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(underestimates_masks_american=1)))
h1b.underestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(underestimates_masks_american=1)),level=.995)
h1b.overestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(overestimates_masks_american=1)))
h1b.overestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(overestimates_masks_american=1)),level=.995)
h2b.underestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(underestimates_masks_copartisans=1)))
h2b.overestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(overestimates_masks_copartisans=1)))
h2b.underestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(underestimates_masks_copartisans=1)),level=.995)
h2b.overestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(overestimates_masks_copartisans=1)),level=.995)
rq1.underestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(underestimates_masks_contrapartisans=1)))
rq1.overestimators.partisans <- summary(margins(Reg.h1b.partisans,at=list(overestimates_masks_contrapartisans=1)))
rq1.underestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(underestimates_masks_contrapartisans=1)),level=.995)
rq1.overestimators.partisans.strict <- confint(margins(Reg.h1b.partisans,at=list(overestimates_masks_contrapartisans=1)),level=.995)
Reg.h1b.partisans.strict <- confint(Reg.h1b.partisans,level=.995)
newdata2 <- data.frame(estimation=factor(c("Underestimated","Accurate","Overestimated","Underestimated","Accurate","Overestimated","Underestimated","Accurate","Overestimated"),levels=c("Overestimated","Accurate","Underestimated")),
                       treatment=factor(c("American","American","American","Co-partisan","Co-partisan","Co-partisan","Out-partisan","Out-partisan","Out-partisan"),levels=c("American","Co-partisan","Out-partisan")),
                       effect=c(h1b.underestimators.partisans$AME[1],Reg.h1b.partisans$coefficients[2],h1b.overestimators.partisans$AME[1],h2b.underestimators.partisans$AME[3],Reg.h1b.partisans$coefficients[5],h2b.overestimators.partisans$AME[3],rq1.underestimators.partisans$AME[2],Reg.h1b.partisans$coefficients[8],rq1.overestimators.partisans$AME[2]),
                       lci=c(h1b.underestimators.partisans$lower[1],Reg.h1b.partisans$conf.low[2],h1b.overestimators.partisans$lower[1],h2b.underestimators.partisans$lower[3],Reg.h1b.partisans$conf.low[5],h2b.overestimators.partisans$lower[3],rq1.underestimators.partisans$lower[2],Reg.h1b.partisans$conf.low[8],rq1.overestimators.partisans$lower[2]),
                       uci=c(h1b.underestimators.partisans$upper[1],Reg.h1b.partisans$conf.high[2],h1b.overestimators.partisans$upper[1],h2b.underestimators.partisans$upper[3],Reg.h1b.partisans$conf.high[5],h2b.overestimators.partisans$upper[3],rq1.underestimators.partisans$upper[2],Reg.h1b.partisans$conf.high[8],rq1.overestimators.partisans$upper[2]),
                       lci_strict=c(h1b.underestimators.partisans.strict[1,1],Reg.h1b.partisans.strict[2,1],h1b.overestimators.partisans.strict[1,1],h2b.underestimators.partisans.strict[4,1],Reg.h1b.partisans.strict[5,1],h2b.overestimators.partisans.strict[4,1],rq1.underestimators.partisans.strict[7,1],Reg.h1b.partisans.strict[8,1],rq1.overestimators.partisans.strict[7,1]),
                       uci_strict=c(h1b.underestimators.partisans.strict[1,2],Reg.h1b.partisans.strict[2,2],h1b.overestimators.partisans.strict[1,2],h2b.underestimators.partisans.strict[4,2],Reg.h1b.partisans.strict[5,2],h2b.overestimators.partisans.strict[4,2],rq1.underestimators.partisans.strict[7,2],Reg.h1b.partisans.strict[8,2],rq1.overestimators.partisans.strict[7,2]))

p.h1b <- ggplot(newdata2,aes(x=estimation,y=effect)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=.8,width=0) + geom_errorbar(aes(ymin=lci_strict,ymax=uci_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + facet_grid(cols=vars(treatment)) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53)) + xlab("Prior norm estimation")

# Figure 3
## Regression
Reg.rq2.partisans <- lm_robust(mask_effective~american_treatment+copartisan_treatment+contrapartisan_treatment+pid3+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2)
## Equivalence Tests
diffmeans.rq2.partisans <- data.partisans %>% group_by(treatment) %>% summarise(mean=mean(mask_effective,na.rm=T),sd=sd(mask_effective,na.rm=T),n=n()) %>% drop_na()
### Control vs. Co-Partisan
TOSTtwo.raw(m1=diffmeans.rq2.partisans$mean[2],m2=diffmeans.rq2.partisans$mean[3],sd1=diffmeans.rq2.partisans$sd[2],sd2=diffmeans.rq2.partisans$sd[3],n1=diffmeans.rq2.partisans$n[2],n2=diffmeans.rq2.partisans$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Control vs. Out-Partisan
TOSTtwo.raw(m1=diffmeans.rq2.partisans$mean[4],m2=diffmeans.rq2.partisans$mean[3],sd1=diffmeans.rq2.partisans$sd[4],sd2=diffmeans.rq2.partisans$sd[3],n1=diffmeans.rq2.partisans$n[4],n2=diffmeans.rq2.partisans$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
## Figure
Reg.rq2.partisans.strict <- confint(Reg.rq2.partisans,level=.995)
newdata3 <- data.frame(treatment=factor(c("American","Copartisan","Contrapartisan"),levels=c("Contrapartisan","Copartisan","American")),
                       effect=c(Reg.rq2.partisans$coefficients[2],Reg.rq2.partisans$coefficients[3],Reg.rq2.partisans$coefficients[4]),
                       lci=c(Reg.rq2.partisans$conf.low[2],Reg.rq2.partisans$conf.low[3],Reg.rq2.partisans$conf.low[4]),
                       uci=c(Reg.rq2.partisans$conf.high[2],Reg.rq2.partisans$conf.high[3],Reg.rq2.partisans$conf.high[4]),
                       lci_strict=c(Reg.rq2.partisans.strict[2,1],Reg.rq2.partisans.strict[3,1],Reg.rq2.partisans.strict[4,1]),
                       uci_strict=c(Reg.rq2.partisans.strict[2,2],Reg.rq2.partisans.strict[3,2],Reg.rq2.partisans.strict[4,2]))

p.rq2 <- ggplot(newdata3,aes(x=treatment,y=effect)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=.8,width=0) + geom_errorbar(aes(ymin=lci_strict,ymax=uci_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53))

# Figure 4
## Regression (Wearing)
Reg.rq3.wear <- lm_robust(wear_mask~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+male+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2)
Reg.rq3.wear.strict <- confint(Reg.rq3.wear,level=.995)
rq3.rep.wear <- summary(margins(Reg.rq3.wear,at=list(republican=1)))
rq3.rep.wear.strict <- confint(margins(Reg.rq3.wear,at=list(republican=1)),level=.995)
## Standardized Effect Size for American Treatent among Republicans
Reg.rq3.wear.gop <- lm_robust(wear_mask~american_treatment+copartisan_treatment+contrapartisan_treatment+male+ideo7+health_trust,data=data,subset=pid3==2)
Reg.rq3.wear.gop$coefficients[2]/sd(subset(data,pid3==2)$wear_mask)
## Figure
### Regression for Effectiveness
Reg.rq3.effective <- lm_robust(mask_effective~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2)
Reg.rq3.effective.strict <- confint(Reg.rq3.effective,level=.995)
rq3.rep.effective <- summary(margins(Reg.rq3.effective,at=list(republican=1)))
rq3.rep.effective.strict <- confint(margins(Reg.rq3.effective,at=list(republican=1)),level=.995)
### Descriptives Generation
rq3.data <- subset(data,pid3==0 | pid3==2)
rq3.wear.data <- rq3.data %>% group_by(pid3,wear_mask) %>% summarise(n=n())
rq3.wear.data$prop <- c(13/1574,10/1574,46/1574,265/1574,1240/1574,59/945,101/945,141/945,261/945,383/945)
rq3.wear.data$Party <- c("Democratic","Democratic","Democratic","Democratic","Democratic","Republican","Republican","Republican","Republican","Republican")
rq3.effective.data <- rq3.data %>% group_by(pid3,mask_effective) %>% summarise(n=n())
rq3.effective.data$prop <- c(28/1574,26/1574,199/1574,1321/1574,111/945,161/945,337/945,336/945)
rq3.effective.data$Party <- c("Democratic","Democratic","Democratic","Democratic","Republican","Republican","Republican","Republican")
rq3.data$Party <- factor(ifelse(rq3.data$pid3==2,"Republican","Democrat"),levels=c("Democrat","Republican"))

p.rq3.1 <- ggplot(rq3.wear.data,aes(x=wear_mask,y=prop,group=Party,fill=Party)) + scale_fill_manual(values=c("Blue","Red")) +  geom_bar(stat="identity",position="dodge") + theme_classic() + ylab("Proportion") + scale_x_continuous(name="Mask-wearing intentions",breaks=c(1,2,3,4,5),labels=c("Not at all","","","","All the time")) + theme(legend.position="none") + ylim(0,1)

newdata4 <- data.frame(treatment=factor(c("American","Copartisan","Contrapartisan","American","Copartisan","Contrapartisan"),levels=c("Contrapartisan","Copartisan","American")),
                       Party=factor(c("Democrat","Democrat","Democrat","Republican","Republican","Republican"),levels=c("Democrat","Republican")),
                       effect_wear=c(Reg.rq3.wear$coefficients[2],Reg.rq3.wear$coefficients[4],Reg.rq3.wear$coefficients[5],rq3.rep.wear$AME[1],rq3.rep.wear$AME[3],rq3.rep.wear$AME[2]),
                       lci_wear=c(Reg.rq3.wear$conf.low[2],Reg.rq3.wear$conf.low[4],Reg.rq3.wear$conf.low[5],rq3.rep.wear$lower[1],rq3.rep.wear$lower[3],rq3.rep.wear$lower[2]),
                       uci_wear=c(Reg.rq3.wear$conf.high[2],Reg.rq3.wear$conf.high[4],Reg.rq3.wear$conf.high[5],rq3.rep.wear$upper[1],rq3.rep.wear$upper[3],rq3.rep.wear$upper[2]),
                       lci_wear_strict=c(Reg.rq3.wear.strict[2,1],Reg.rq3.wear.strict[4,1],Reg.rq3.wear.strict[5,1],rq3.rep.wear.strict[1,1],rq3.rep.wear.strict[3,1],rq3.rep.wear.strict[4,1]),
                       uci_wear_strict=c(Reg.rq3.wear.strict[2,2],Reg.rq3.wear.strict[4,2],Reg.rq3.wear.strict[5,2],rq3.rep.wear.strict[1,2],rq3.rep.wear.strict[3,2],rq3.rep.wear.strict[4,2]),
                       effect_effective=c(Reg.rq3.effective$coefficients[2],Reg.rq3.effective$coefficients[4],Reg.rq3.effective$coefficients[5],rq3.rep.effective$AME[1],rq3.rep.effective$AME[3],rq3.rep.effective$AME[2]),
                       lci_effective=c(Reg.rq3.effective$conf.low[2],Reg.rq3.effective$conf.low[4],Reg.rq3.effective$conf.low[5],rq3.rep.effective$lower[1],rq3.rep.effective$lower[3],rq3.rep.effective$lower[2]),
                       uci_effective=c(Reg.rq3.effective$conf.high[2],Reg.rq3.effective$conf.high[4],Reg.rq3.effective$conf.high[5],rq3.rep.effective$upper[1],rq3.rep.effective$upper[3],rq3.rep.effective$upper[2]),
                       lci_effective_strict=c(Reg.rq3.effective.strict[2,1],Reg.rq3.effective.strict[4,1],Reg.rq3.effective.strict[5,1],rq3.rep.effective.strict[1,1],rq3.rep.effective.strict[3,1],rq3.rep.effective.strict[4,1]),
                       uci_effective_strict=c(Reg.rq3.effective.strict[2,2],Reg.rq3.effective.strict[4,2],Reg.rq3.effective.strict[5,2],rq3.rep.effective.strict[1,2],rq3.rep.effective.strict[3,2],rq3.rep.effective.strict[4,2]))

p.rq3.3 <- ggplot(newdata4,aes(x=treatment,y=effect_wear,group=Party,color=Party)) + scale_color_manual(values=c("Blue","Red")) + coord_flip() + geom_point(position=position_dodge(0.25)) + geom_errorbar(aes(ymin=lci_wear,ymax=uci_wear),size=.8,width=0,position=position_dodge(0.25)) + geom_errorbar(aes(ymin=lci_wear_strict,ymax=uci_wear_strict),size=.4,width=0,position=position_dodge(0.25)) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Black") + theme(legend.position="none") + scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53))

p3.1 <- ggarrange(p.rq3.1,p.rq3.3,nrow=1,common.legend=T,legend="bottom")

# Table A1
a1.m1 <- lm_robust(wear_mask~american_treatment+copartisan_treatment+contrapartisan_treatment,data=data,subset=pid3==0 | pid3==2)
a1.m2 <- lm_robust(wear_mask~american_treatment+copartisan_treatment+contrapartisan_treatment+male+pid3+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2)

# Table A2
a2.m1 <- lm_robust(wear_mask~american_treatment,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)
a2.m2 <- lm_robust(wear_mask~american_treatment+male+pid3+ideo7+health_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)
## Standardized Effect Size for American Treatment
a2.m2$coefficients[2]/sd(subset(data,copartisan_treatment==0 & contrapartisan_treatment==0)$wear_mask)

# Table A3
a3.m1 <- lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans),data=data)
a3.m2 <- lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans)+male+pid3+ideo7+health_trust,data=data)

# Table A4
a4.m1 <- lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american),data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)
a4.m2 <- lm_robust(wear_mask~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+male+pid3+ideo7+health_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)

# Table A5
a5.m1 <- lm_robust(mask_effective~american_treatment+copartisan_treatment+contrapartisan_treatment,data=data,subset=pid3==0 | pid3==2)
a5.m2 <- lm_robust(mask_effective~american_treatment+copartisan_treatment+contrapartisan_treatment+pid3+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2)

# Table A6
a6.m1 <- lm_robust(mask_effective~american_treatment,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)
a6.m2 <- lm_robust(mask_effective~american_treatment+pid3+ideo7+health_trust+media_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0)
## Equivalence Test for American vs. Control
data.full <- subset(data,copartisan_treatment==0 & contrapartisan_treatment==0)
diffmeans.rq2.full <- data.full %>% group_by(treatment) %>% summarise(mean=mean(mask_effective,na.rm=T),sd=sd(mask_effective,na.rm=T),n=n()) %>% drop_na()
TOSTtwo.raw(m1=diffmeans.rq2.full$mean[1],m2=diffmeans.rq2.full$mean[2],sd1=diffmeans.rq2.full$sd[1],sd2=diffmeans.rq2.full$sd[2],n1=diffmeans.rq2.full$n[1],n2=diffmeans.rq2.full$n[2],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 

# Table A7
a7.m1 <- lm_robust(wear_mask~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican),data=data,subset=pid3==0 | pid3==2)
a7.m2 <- lm_robust(wear_mask~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+male+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2)
a7.m3 <- lm_robust(mask_effective~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican),data=data,subset=pid3==0 | pid3==2)
a7.m4 <- lm_robust(mask_effective~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2)

# Figure A1 (See code for Figure 3 for how it was created)
p.rq3.2 <- ggplot(rq3.effective.data,aes(x=mask_effective,y=prop,group=Party,fill=Party)) + scale_fill_manual(values=c("Blue","Red")) +  geom_bar(stat="identity",position="dodge") + theme_classic() + ylab("Proportion") + scale_x_continuous(name="Mask effectiveness (accuracy)",breaks=c(1,2,3,4),labels=c("Not at all","","","Very")) + theme(legend.position="none") + ylim(0,1)
p.rq3.4 <- ggplot(newdata4,aes(x=treatment,y=effect_effective,group=Party,color=Party)) + scale_color_manual(values=c("Blue","Red")) + coord_flip() + geom_point(position=position_dodge(0.25)) + geom_errorbar(aes(ymin=lci_effective,ymax=uci_effective),size=.8,width=0,position=position_dodge(0.25)) + geom_errorbar(aes(ymin=lci_effective_strict,ymax=uci_effective_strict),size=.4,width=0,position=position_dodge(0.25)) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Black") + theme(legend.position="none") + scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53))
p3.2 <- ggarrange(p.rq3.2,p.rq3.4,nrow=1,common.legend=T,legend="bottom")

# Table A8
## Regressions
a8.m1 <- lm_robust(affective_polarization~american_treatment+copartisan_treatment+contrapartisan_treatment,data=data)
a8.m2 <- lm_robust(affective_polarization~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican),data=data)
a8.m3 <- lm_robust(affective_polarization~american_treatment+copartisan_treatment+contrapartisan_treatment+ideo7+polinterest,data=data)
a8.m4 <- lm_robust(affective_polarization~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+polinterest,data=data)
## Equivalence Tests
diffmeans.rq5.1 <- data.partisans %>% group_by(treatment) %>% summarise(mean=mean(affective_polarization,na.rm=T),
                                                                        sd=sd(affective_polarization,na.rm=T),
                                                                        n=n()) %>% drop_na()

diffmeans.rq5.2 <- data.partisans %>% group_by(treatment,pid3) %>% summarise(mean=mean(affective_polarization,na.rm=T),
                                                                             sd=sd(affective_polarization,na.rm=T),
                                                                             n=n()) %>% drop_na()

### American vs. Control
TOSTtwo.raw(m1=diffmeans.rq5.1$mean[1],m2=diffmeans.rq5.1$mean[3],sd1=diffmeans.rq5.1$sd[1],sd2=diffmeans.rq5.1$sd[3],n1=diffmeans.rq5.1$n[1],n2=diffmeans.rq5.1$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Co-Partisan vs. Control
TOSTtwo.raw(m1=diffmeans.rq5.1$mean[2],m2=diffmeans.rq5.1$mean[3],sd1=diffmeans.rq5.1$sd[2],sd2=diffmeans.rq5.1$sd[3],n1=diffmeans.rq5.1$n[2],n2=diffmeans.rq5.1$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Out-Partisan vs. Control
TOSTtwo.raw(m1=diffmeans.rq5.1$mean[4],m2=diffmeans.rq5.1$mean[3],sd1=diffmeans.rq5.1$sd[4],sd2=diffmeans.rq5.1$sd[3],n1=diffmeans.rq5.1$n[4],n2=diffmeans.rq5.1$n[3],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
## American vs. Control, Democrats
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[1],m2=diffmeans.rq5.2$mean[5],sd1=diffmeans.rq5.2$sd[1],sd2=diffmeans.rq5.2$sd[5],n1=diffmeans.rq5.2$n[1],n2=diffmeans.rq5.2$n[5],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Co-Partisan vs. Control, Democrats
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[3],m2=diffmeans.rq5.2$mean[5],sd1=diffmeans.rq5.2$sd[3],sd2=diffmeans.rq5.2$sd[5],n1=diffmeans.rq5.2$n[3],n2=diffmeans.rq5.2$n[5],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Out-Partisan vs. Control, Democrats
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[7],m2=diffmeans.rq5.2$mean[5],sd1=diffmeans.rq5.2$sd[7],sd2=diffmeans.rq5.2$sd[5],n1=diffmeans.rq5.2$n[7],n2=diffmeans.rq5.2$n[5],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### American vs. Control, Republicans
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[2],m2=diffmeans.rq5.2$mean[6],sd1=diffmeans.rq5.2$sd[2],sd2=diffmeans.rq5.2$sd[6],n1=diffmeans.rq5.2$n[2],n2=diffmeans.rq5.2$n[6],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Co-Partisan vs. Control, Republicans
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[4],m2=diffmeans.rq5.2$mean[6],sd1=diffmeans.rq5.2$sd[4],sd2=diffmeans.rq5.2$sd[6],n1=diffmeans.rq5.2$n[4],n2=diffmeans.rq5.2$n[6],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 
### Out-Partisan vs. Control, Republicans
TOSTtwo.raw(m1=diffmeans.rq5.2$mean[8],m2=diffmeans.rq5.2$mean[6],sd1=diffmeans.rq5.2$sd[8],sd2=diffmeans.rq5.2$sd[6],n1=diffmeans.rq5.2$n[8],n2=diffmeans.rq5.2$n[6],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F) 

# Figure A2 (based on Table A8)
Reg.rq5.basic <- lm_robust(affective_polarization~american_treatment+copartisan_treatment+contrapartisan_treatment+ideo7+polinterest,data=data)
Reg.rq5.basic.strict <- confint(Reg.rq5.basic,level=.995)
newdata6 <- data.frame(treatment=factor(c("American","Copartisan","Contrapartisan"),levels=c("Contrapartisan","Copartisan","American")),
                       effect=c(Reg.rq5.basic$coefficients[2],Reg.rq5.basic$coefficients[3],Reg.rq5.basic$coefficients[4]),
                       lci=c(Reg.rq5.basic$conf.low[2],Reg.rq5.basic$conf.low[3],Reg.rq5.basic$conf.low[4]),
                       uci=c(Reg.rq5.basic$conf.high[2],Reg.rq5.basic$conf.high[3],Reg.rq5.basic$conf.high[4]),
                       lci_strict=c(Reg.rq5.basic.strict[2,1],Reg.rq5.basic.strict[3,1],Reg.rq5.basic.strict[4,1]),
                       uci_strict=c(Reg.rq5.basic.strict[2,2],Reg.rq5.basic.strict[3,2],Reg.rq5.basic.strict[4,2]))

p.rq5 <- ggplot(newdata6,aes(x=treatment,y=effect)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci,ymax=uci),size=.8,width=0) + geom_errorbar(aes(ymin=lci_strict,ymax=uci_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + ylab("Treatment effect")

# Table A9
## Regressions
a9.m1 <- lm_robust(wear_mask~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment),data=data.partisans)
a9.m2 <- lm_robust(wear_mask~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+male+pid3+ideo7+health_trust,data=data.partisans)
a9.m3 <- lm_robust(mask_effective~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment),data=data.partisans)
a9.m4 <- lm_robust(mask_effective~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+pid3+ideo7+health_trust+media_trust,data=data.partisans)
## Interaction F-tests
### Wearing
linearHypothesis(a9.m2,c("factcheck_treatment=0","W3factcheck_treatment=0","american_treatment:factcheck_treatment=0","american_treatment:W3factcheck_treatment=0","factcheck_treatment:copartisan_treatment=0","W3factcheck_treatment:copartisan_treatment=0","factcheck_treatment:contrapartisan_treatment=0","W3factcheck_treatment:contrapartisan_treatment=0"),test="F") ## p=.2782, fail to reject null
### Effectiveness
linearHypothesis(a9.m4,c("factcheck_treatment=0","W3factcheck_treatment=0","american_treatment:factcheck_treatment=0","american_treatment:W3factcheck_treatment=0","factcheck_treatment:copartisan_treatment=0","W3factcheck_treatment:copartisan_treatment=0","factcheck_treatment:contrapartisan_treatment=0","W3factcheck_treatment:contrapartisan_treatment=0"),test="F") ## p=.014, interpret coefficients
## Equivalence Tests
diffmeans.rq4.2 <- data.partisans %>% group_by(treatment,factcheck_treatment,W3factcheck_treatment) %>% summarise(mean=mean(mask_effective,na.rm=T),sd=sd(mask_effective,na.rm=T),n=n()) %>% drop_na()
diffmeans.rq4.4 <- data.partisans %>% group_by(treatment,factcheck_treatment) %>% summarise(mean=mean(mask_effective,na.rm=T),sd=sd(mask_effective,na.rm=T),n=n()) %>% drop_na()
diffmeans.rq4.6 <- data.partisans %>% group_by(treatment,W3factcheck_treatment) %>% summarise(mean=mean(mask_effective,na.rm=T),sd=sd(mask_effective,na.rm=T),n=n()) %>% drop_na()
### Control vs. American, No Fact Checks
TOSTtwo.raw(m1=diffmeans.rq4.2$mean[5],m2=diffmeans.rq4.2$mean[9],sd1=diffmeans.rq4.2$sd[5],sd2=diffmeans.rq4.2$sd[9],n1=diffmeans.rq4.2$n[5],n2=diffmeans.rq4.2$n[9],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F)
### Control vs. American, W3 Fact-Check
TOSTtwo.raw(m1=diffmeans.rq4.4$mean[4],m2=diffmeans.rq4.4$mean[6],sd1=diffmeans.rq4.4$sd[4],sd2=diffmeans.rq4.4$sd[6],n1=diffmeans.rq4.4$n[4],n2=diffmeans.rq4.4$n[6],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F)
### Control vs. American, Both Fact-Checks
TOSTtwo.raw(m1=diffmeans.rq4.6$mean[4],m2=diffmeans.rq4.6$mean[6],sd1=diffmeans.rq4.6$sd[4],sd2=diffmeans.rq4.6$sd[6],n1=diffmeans.rq4.6$n[4],n2=diffmeans.rq4.6$n[6],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F)
### Control vs. Out-Partisan, W2 Fact-Check
TOSTtwo.raw(m1=diffmeans.rq4.2$mean[13],m2=diffmeans.rq4.2$mean[9],sd1=diffmeans.rq4.2$sd[13],sd2=diffmeans.rq4.2$sd[9],n1=diffmeans.rq4.2$n[13],n2=diffmeans.rq4.2$n[9],low_eqbound=-0.1,high_eqbound=0.1,alpha=0.05,var.equal=F)

# Figure A3
## Regressions
data.partisans$w2only <- ifelse(data.partisans$factcheck_treatment==1 & data.partisans$W3factcheck_treatment==0,1,0)
data.partisans$w3only <- ifelse(data.partisans$factcheck_treatment==0 & data.partisans$W3factcheck_treatment==1,1,0)
data.partisans$both <- ifelse(data.partisans$factcheck_treatment==1 & data.partisans$W3factcheck_treatment==1,1,0)
Reg.rq4.wear <- lm_robust(wear_mask~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+male+pid3+ideo7+health_trust,data=data.partisans)
Reg.rq4.wear.strict <- confint(Reg.rq4.wear,level=.995)
Reg.rq4.effective <- lm_robust(mask_effective~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+pid3+ideo7+health_trust+media_trust,data=data.partisans)
Reg.rq4.effective.strict <- confint(Reg.rq4.effective,level=.995)
Reg.rq4.wear.allinteractions <- lm_robust(wear_mask~(american_treatment*w2only)+(american_treatment*w3only)+(american_treatment*both)+(copartisan_treatment*w2only)+(copartisan_treatment*w3only)+(copartisan_treatment*both)+(contrapartisan_treatment*w2only)+(contrapartisan_treatment*w3only)+(contrapartisan_treatment*both),data=data.partisans)
Reg.rq4.wear.allinteractions.strict <- confint(Reg.rq4.wear.allinteractions,level=.995)
Reg.rq4.effective.allinteractions <- lm_robust(mask_effective~(american_treatment*w2only)+(american_treatment*w3only)+(american_treatment*both)+(copartisan_treatment*w2only)+(copartisan_treatment*w3only)+(copartisan_treatment*both)+(contrapartisan_treatment*w2only)+(contrapartisan_treatment*w3only)+(contrapartisan_treatment*both)+pid3+ideo7+health_trust+media_trust,data=data.partisans)
Reg.rq4.effective.allinteractions.strict <- confint(Reg.rq4.effective.allinteractions,level=.995)
rq4.wear.nofactchecks <- summary(margins(Reg.rq4.wear,at=list(factcheck_treatment=0,W3factcheck_treatment=0)))
rq4.wear.w2factcheck <- summary(margins(Reg.rq4.wear,at=list(factcheck_treatment=1)))
rq4.wear.w3factcheck <- summary(margins(Reg.rq4.wear,at=list(W3factcheck_treatment=1)))
rq4.wear.bothfactcheck <- summary(margins(Reg.rq4.wear,at=list(factcheck_treatment=1,W3factcheck_treatment=1)))
rq4.wear.nofactchecks.strict <- confint(margins(Reg.rq4.wear,at=list(factcheck_treatment=0,W3factcheck_treatment=0)),level=.995)
rq4.wear.w2factcheck.strict <- confint(margins(Reg.rq4.wear,at=list(factcheck_treatment=1)),level=.995)
rq4.wear.w3factcheck.strict <- confint(margins(Reg.rq4.wear,at=list(W3factcheck_treatment=1)),level=.995)
rq4.wear.bothfactcheck.strict <- confint(margins(Reg.rq4.wear,at=list(factcheck_treatment=1,W3factcheck_treatment=1)),level=.995)
rq4.effective.nofactchecks <- summary(margins(Reg.rq4.effective,at=list(factcheck_treatment=0,W3factcheck_treatment=0)))
rq4.effective.w2factcheck <- summary(margins(Reg.rq4.effective,at=list(factcheck_treatment=1)))
rq4.effective.w3factcheck <- summary(margins(Reg.rq4.effective,at=list(W3factcheck_treatment=1)))
rq4.effective.bothfactcheck <- summary(margins(Reg.rq4.effective,at=list(factcheck_treatment=1,W3factcheck_treatment=1)))
rq4.effective.nofactchecks.strict <- confint(margins(Reg.rq4.effective,at=list(factcheck_treatment=0,W3factcheck_treatment=0)),level=.995)
rq4.effective.w2factcheck.strict <- confint(margins(Reg.rq4.effective,at=list(factcheck_treatment=1)),level=.995)
rq4.effective.w3factcheck.strict <- confint(margins(Reg.rq4.effective,at=list(W3factcheck_treatment=1)),level=.995)
rq4.effective.bothfactcheck.strict <- confint(margins(Reg.rq4.effective,at=list(factcheck_treatment=1,W3factcheck_treatment=1)),level=.995)
## Figure
newdata5 <- data.frame(treatment=factor(c("American","Copartisan","Contrapartisan","American","Copartisan","Contrapartisan","American","Copartisan","Contrapartisan","American","Copartisan","Contrapartisan"),levels=c("Contrapartisan","Copartisan","American")),
                       factcheck=factor(c("None","None","None","W2 Only","W2 Only","W2 Only","W3 Only","W3 Only","W3 Only","W2 and W3","W2 and W3","W2 and W3"),levels=c("None","W2 Only","W3 Only","W2 and W3")),
                       effect_wear=c(Reg.rq4.wear$coefficients[2],Reg.rq4.wear$coefficients[5],Reg.rq4.wear$coefficients[6],rq4.wear.w2factcheck$AME[1],rq4.wear.w2factcheck$AME[3],rq4.wear.w2factcheck$AME[2],rq4.wear.w3factcheck$AME[1],rq4.wear.w3factcheck$AME[3],rq4.wear.w3factcheck$AME[2],rq4.wear.bothfactcheck$AME[1],rq4.wear.bothfactcheck$AME[3],rq4.wear.bothfactcheck$AME[2]),
                       lci_wear=c(Reg.rq4.wear$conf.low[2],Reg.rq4.wear$conf.low[5],Reg.rq4.wear$conf.low[6],rq4.wear.w2factcheck$lower[1],rq4.wear.w2factcheck$lower[3],rq4.wear.w2factcheck$lower[2],rq4.wear.w3factcheck$lower[1],rq4.wear.w3factcheck$lower[3],rq4.wear.w3factcheck$lower[2],rq4.wear.bothfactcheck$lower[1],rq4.wear.bothfactcheck$lower[3],rq4.wear.bothfactcheck$lower[2]),
                       uci_wear=c(Reg.rq4.wear$conf.high[2],Reg.rq4.wear$conf.high[5],Reg.rq4.wear$conf.high[6],rq4.wear.w2factcheck$upper[1],rq4.wear.w2factcheck$upper[3],rq4.wear.w2factcheck$upper[2],rq4.wear.w3factcheck$upper[1],rq4.wear.w3factcheck$upper[3],rq4.wear.w3factcheck$upper[2],rq4.wear.bothfactcheck$upper[1],rq4.wear.bothfactcheck$upper[3],rq4.wear.bothfactcheck$upper[2]),
                       lci_wear_strict=c(Reg.rq4.wear.strict[2,1],Reg.rq4.wear.strict[5,1],Reg.rq4.wear.strict[6,1],rq4.wear.w2factcheck.strict[1,1],rq4.wear.w2factcheck.strict[4,1],rq4.wear.w2factcheck.strict[5,1],rq4.wear.w3factcheck.strict[1,1],rq4.wear.w3factcheck.strict[4,1],rq4.wear.w3factcheck.strict[5,1],rq4.wear.bothfactcheck.strict[1,1],rq4.wear.bothfactcheck.strict[4,1],rq4.wear.bothfactcheck.strict[5,1]),
                       uci_wear_strict=c(Reg.rq4.wear.strict[2,2],Reg.rq4.wear.strict[5,2],Reg.rq4.wear.strict[6,2],rq4.wear.w2factcheck.strict[1,2],rq4.wear.w2factcheck.strict[4,2],rq4.wear.w2factcheck.strict[5,2],rq4.wear.w3factcheck.strict[1,2],rq4.wear.w3factcheck.strict[4,2],rq4.wear.w3factcheck.strict[5,2],rq4.wear.bothfactcheck.strict[1,2],rq4.wear.bothfactcheck.strict[4,2],rq4.wear.bothfactcheck.strict[5,2]),
                       effect_effective=c(Reg.rq4.effective$coefficients[2],Reg.rq4.effective$coefficients[5],Reg.rq4.effective$coefficients[6],rq4.effective.w2factcheck$AME[1],rq4.effective.w2factcheck$AME[3],rq4.effective.w2factcheck$AME[2],rq4.effective.w3factcheck$AME[1],rq4.effective.w3factcheck$AME[3],rq4.effective.w3factcheck$AME[2],rq4.effective.bothfactcheck$AME[1],rq4.effective.bothfactcheck$AME[3],rq4.effective.bothfactcheck$AME[2]),
                       lci_effective=c(Reg.rq4.effective$conf.low[2],Reg.rq4.effective$conf.low[5],Reg.rq4.effective$conf.low[6],rq4.effective.w2factcheck$lower[1],rq4.effective.w2factcheck$lower[3],rq4.effective.w2factcheck$lower[2],rq4.effective.w3factcheck$lower[1],rq4.effective.w3factcheck$lower[3],rq4.effective.w3factcheck$lower[2],rq4.effective.bothfactcheck$lower[1],rq4.effective.bothfactcheck$lower[3],rq4.effective.bothfactcheck$lower[2]),
                       uci_effective=c(Reg.rq4.effective$conf.high[2],Reg.rq4.effective$conf.high[5],Reg.rq4.effective$conf.high[6],rq4.effective.w2factcheck$upper[1],rq4.effective.w2factcheck$upper[3],rq4.effective.w2factcheck$upper[2],rq4.effective.w3factcheck$upper[1],rq4.effective.w3factcheck$upper[3],rq4.effective.w3factcheck$upper[2],rq4.effective.bothfactcheck$upper[1],rq4.effective.bothfactcheck$upper[3],rq4.effective.bothfactcheck$upper[2]),
                       lci_effective_strict=c(Reg.rq4.effective.strict[2,1],Reg.rq4.effective.strict[5,1],Reg.rq4.effective.strict[6,1],rq4.effective.w2factcheck.strict[1,1],rq4.effective.w2factcheck.strict[4,1],rq4.effective.w2factcheck.strict[5,1],rq4.effective.w3factcheck.strict[1,1],rq4.effective.w3factcheck.strict[4,1],rq4.effective.w3factcheck.strict[5,1],rq4.effective.bothfactcheck.strict[1,1],rq4.effective.bothfactcheck.strict[4,1],rq4.effective.bothfactcheck.strict[5,1]),
                       uci_effective_strict=c(Reg.rq4.effective.strict[2,2],Reg.rq4.effective.strict[5,2],Reg.rq4.effective.strict[6,2],rq4.effective.w2factcheck.strict[1,2],rq4.effective.w2factcheck.strict[4,2],rq4.effective.w2factcheck.strict[5,2],rq4.effective.w3factcheck.strict[1,2],rq4.effective.w3factcheck.strict[4,2],rq4.effective.w3factcheck.strict[5,2],rq4.effective.bothfactcheck.strict[1,2],rq4.effective.bothfactcheck.strict[4,2],rq4.effective.bothfactcheck.strict[5,2]))

p.rq4.1 <- ggplot(newdata5,aes(x=treatment,y=effect_wear)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci_wear,ymax=uci_wear),size=.8,width=0) + geom_errorbar(aes(ymin=lci_wear_strict,ymax=uci_wear_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + facet_grid(rows=vars(factcheck))+ scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53)) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Mask-wearing intentions")
p.rq4.2 <- ggplot(newdata5,aes(x=treatment,y=effect_effective)) + coord_flip() + geom_point() + geom_errorbar(aes(ymin=lci_effective,ymax=uci_effective),size=.8,width=0) + geom_errorbar(aes(ymin=lci_effective_strict,ymax=uci_effective_strict),size=.4,width=0) + theme_classic() + geom_hline(yintercept=0,linetype="dashed",color="Red") + facet_grid(rows=vars(factcheck))+ scale_x_discrete(name="Norm treatment",breaks=c("Contrapartisan","Copartisan","American"),labels=c("Out-partisan","Co-partisan","American")) + scale_y_continuous(name="Treatment effect",limits=c(-0.29,0.53)) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("Mask effectiveness")
p4.1 <- grid.arrange(p.rq4.1,p.rq4.2,nrow=1)

# Table A10
a10.m1 <- polr(as.factor(wear_mask)~american_treatment+copartisan_treatment+contrapartisan_treatment,data=data,subset=pid3==0 | pid3==2,Hess=T)
a10.m2 <- polr(as.factor(wear_mask)~american_treatment+copartisan_treatment+contrapartisan_treatment+male+pid3+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2,Hess=T)

# Table A11
## Set Up Synthetic data
newdata.h1a.partisans <- data.frame(american_treatment=c(0,1,0,0),copartisan_treatment=c(0,0,1,0),contrapartisan_treatment=c(0,0,0,1),male=0,pid3=0.7503,ideo7=3.84,health_trust=2.275)
## Predicted Probabilities
Reg.h1a.partisans.probs <- predict(a10.m2,newdata.h1a.partisans,type="probs")

# Table A12
a12.m1 <- polr(as.factor(wear_mask)~american_treatment,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)
a12.m2 <- polr(as.factor(wear_mask)~american_treatment+male+pid3+ideo7+health_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)

# Table A13
a13.m1 <- polr(as.factor(wear_mask)~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans),data=data,Hess=T)
a13.m2 <- polr(as.factor(wear_mask)~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+(copartisan_treatment*underestimates_masks_copartisans)+(copartisan_treatment*overestimates_masks_copartisans)+(contrapartisan_treatment*underestimates_masks_contrapartisans)+(contrapartisan_treatment*overestimates_masks_contrapartisans)+male+pid3+ideo7+health_trust,data=data,Hess=T)

# Table A14
a14.m1 <- polr(as.factor(wear_mask)~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american),data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)
a14.m2 <- polr(as.factor(wear_mask)~(american_treatment*underestimates_masks_american)+(american_treatment*overestimates_masks_american)+male+pid3+ideo7+health_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)

# Table A15
## Partisans
### Prepare Synthetic Data
newdata.h1b.partisans <- data.frame(american_treatment=c(0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0),copartisan_treatment=c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0),contrapartisan_treatment=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1),underestimates_masks_american=c(0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0),overestimates_masks_american=c(0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0),underestimates_masks_copartisans=c(0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0),overestimates_masks_copartisans=c(0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0),underestimates_masks_contrapartisans=c(0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0),overestimates_masks_contrapartisans=c(0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1),male=0,pid3=0.750,ideo7=3.84,health_trust=2.275)
### Predict Probabilities
Reg.h1b.partisans.probs <- predict(a13.m2,newdata.h1b.partisans,type="probs")
## Full Sample
### Prepare Synthetic Data
newdata.h1b.full <- data.frame(american_treatment=c(0,0,0,1,1,1),underestimates_masks_american=c(0,1,0,0,1,0),overestimates_masks_american=c(0,0,1,0,0,1),male=0,pid3=0.8104,ideo7=3.919,health_trust=2.222)
### Predict Probabilities
Reg.h1b.full.probs <- predict(a14.m2,newdata.h1b.full,type="probs")

# Table A16
a16.m1 <- polr(as.factor(mask_effective)~american_treatment+copartisan_treatment+contrapartisan_treatment,data=data,subset=pid3==0 | pid3==2,Hess=T)
a16.m2 <- polr(as.factor(mask_effective)~american_treatment+copartisan_treatment+contrapartisan_treatment+pid3+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2,Hess=T)

# Table A17
a17.m1 <- polr(as.factor(mask_effective)~american_treatment,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)
a17.m2 <- polr(as.factor(mask_effective)~american_treatment+pid3+ideo7+health_trust+media_trust,data=data,subset=copartisan_treatment==0 & contrapartisan_treatment==0,Hess=T)

# Table A18
## Partisans
### Prepare Synthetic Data
newdata.rq2.partisans <- data.frame(american_treatment=c(0,1,0,0),copartisan_treatment=c(0,0,1,0),contrapartisan_treatment=c(0,0,0,1),media_trust=1.949,pid3=0.7503,ideo7=3.84,health_trust=2.275)
### Predict Probabilities
Reg.rq2.partisans.probs <- predict(a16.m2,newdata.rq2.partisans,type="probs")
## Full Sample
### Prepare Synthetic Data
newdata.rq2.full <- data.frame(american_treatment=c(0,1),media_trust=1.865,pid3=0.8104,ideo7=3.919,health_trust=2.222)
### Predict Probabilities
Reg.rq2.full.probs <- predict(a17.m2,newdata.rq2.full,type="probs")

# Table A19
a19.m1 <- polr(as.factor(wear_mask)~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican),data=data,subset=pid3==0 | pid3==2,Hess=T)
a19.m2 <- polr(as.factor(wear_mask)~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+male+ideo7+health_trust,data=data,subset=pid3==0 | pid3==2,Hess=T)
a19.m3 <- polr(as.factor(mask_effective)~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican),data=data,subset=pid3==0 | pid3==2,Hess=T)
a19.m4 <- polr(as.factor(mask_effective)~(american_treatment*republican)+(copartisan_treatment*republican)+(contrapartisan_treatment*republican)+ideo7+health_trust+media_trust,data=data,subset=pid3==0 | pid3==2,Hess=T)

# Table A20
## Prepare Synthetic Data
newdata.rq3.wear <- data.frame(american_treatment=c(0,0,1,1,0,0,0,0),copartisan_treatment=c(0,0,0,0,1,1,0,0),republican=c(0,1,0,1,0,1,0,1),contrapartisan_treatment=c(0,0,0,0,0,0,1,1),male=0,pid3=0.7503,ideo7=3.84,health_trust=2.275)
## Predict Probabilities
Reg.rq3.wear.probs <- predict(a19.m2,newdata.rq3.wear,type="probs")

# Table A21
## Prepare Synthetic Data
newdata.rq3.effective <- data.frame(american_treatment=c(0,0,1,1,0,0,0,0),copartisan_treatment=c(0,0,0,0,1,1,0,0),republican=c(0,1,0,1,0,1,0,1),contrapartisan_treatment=c(0,0,0,0,0,0,1,1),media_trust=1.949,pid3=0.7503,ideo7=3.84,health_trust=2.275)
## Predict Probabilities
Reg.rq3.effective.probs <- predict(a19.m4,newdata.rq3.effective,type="probs")

# Table A22
a22.m1 <- polr(as.factor(wear_mask)~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment),data=data.partisans,Hess=T)
a22.m2 <- polr(as.factor(wear_mask)~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+male+pid3+ideo7+health_trust,data=data.partisans,Hess=T)
a22.m3 <- polr(as.factor(mask_effective)~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment),data=data.partisans,Hess=T)
a22.m4 <- polr(as.factor(mask_effective)~(american_treatment*factcheck_treatment)+(american_treatment*W3factcheck_treatment)+(copartisan_treatment*factcheck_treatment)+(copartisan_treatment*W3factcheck_treatment)+(contrapartisan_treatment*factcheck_treatment)+(contrapartisan_treatment*W3factcheck_treatment)+pid3+ideo7+health_trust+media_trust,data=data.partisans,Hess=T)

# Table A23
## Prepare Synthetic Data
newdata.rq4.wear <- data.frame(american_treatment=c(0,0,0,1,1,1,0,0,0,0,0,0),copartisan_treatment=c(0,0,0,0,0,0,1,1,1,0,0,0),contrapartisan_treatment=c(0,0,0,0,0,0,0,0,0,1,1,1),factcheck_treatment=c(0,1,0,0,1,0,0,1,0,0,1,0),W3factcheck_treatment=c(0,0,1,0,0,1,0,0,1,0,0,1),male=0,pid3=0.7503,ideo7=3.84,health_trust=2.275)
## Predict Probabilities
Reg.rq4.wear.probs <- predict(Reg.rq4.wear.ologit,newdata.rq4.wear,type="probs")

# Table A24
## Prepare Synthetic Data
newdata.rq4.effective <- data.frame(american_treatment=c(0,0,0,1,1,1,0,0,0,0,0,0),copartisan_treatment=c(0,0,0,0,0,0,1,1,1,0,0,0),contrapartisan_treatment=c(0,0,0,0,0,0,0,0,0,1,1,1),factcheck_treatment=c(0,1,0,0,1,0,0,1,0,0,1,0),W3factcheck_treatment=c(0,0,1,0,0,1,0,0,1,0,0,1),media_trust=1.949,pid3=0.7503,ideo7=3.84,health_trust=2.275)
## Predict Probabilities
Reg.rq4.effective.probs <- predict(Reg.rq4.effective.ologit,newdata.rq4.effective,type="probs")
